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Abstract. With the emergence in the next few years of a new breed of high power 
laser facilities, it is becoming increasingly important to understand how interacting 
with intense laser pulses affects the bulk properties of a relativistic electron beam. A 
detailed analysis of the radiative cooling of electrons indicates that, classically, equal 
contributions to the phase space contraction occur in the transverse and longitudinal 
directions. In the weakly quantum regime, in addition to an overall reduction in beam 
cooling, this symmetry is broken, leading to significantly less cooling in the longitudinal 
than the transverse directions. By introducing an efficient new technique for studying 
the evolution of a particle distribution, we demonstrate the quantum reduction in beam 
cooling, and find that it depends on the distribution of energy in the laser pulse, rather 
than just the total energy as in the classical case. 
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1 . Introduction 


The emergence over the next few years of a new generation of ultra-high power laser 
facilities, spearheaded by the Extreme Light Infrastructure (ELI) [1], represents a major 
advance in the possibilities afforded by laser technology. In addition to important 
practical applications, these facilities will, for the hrst time, allow investigation of 
qualitatively new physical regimes. Among the hrst effects to be explored will be 
radiation reaction. 

Radiation reaction—the recoil force on an electron due to its emission of radiation— 
remains a contentious area of physics after more than a century of investigation. The 
standard equation describing radiation reaction (the so-called LAD equation, after its 
progenitors Lorentz, Abraham, and Dirac [2, 3, 4]) for a particle of mass m and charge 
q in an electromagnetic held F reads 

= ^ + tA\x’’ = -^F\x^ + rCf “ - xbi^x^), (1) 

m m 

where = —qF'^bX^ is the Lorentz force. Here, the constant r := g^/bvrm is 
the ‘characteristic time’ of the particle| and an overdot denotes diherentiation with 
respect to proper time. Indices are raised and lowered with the metric tensor rj = 
diag(—1,1,1,1), and repeated indices are summed from 0 to 3. The i-orthogonal 
projection := 6^ + x°'Xh ensures that x is orthogonal to x, preserving the 

normalisation condition x°’Xa = — 1 (equivalently the mass shell condition, p“pa = 
where = mx°‘ = (ym, p)). We work in Heaviside-Lorentz units with c = 1. 

Equation (1) may be unpacked and expressed in terms of the 3-momentum as 


% = JE+-^xB|+r7 


dt 


'jm 


i-( 7 Tj+p(^j' 

j dt y dt j ^ \ dt j 


JL ( 

m? \dt J ’ 


( 2 ) 


where 7 = -I- p'^/nF and dj/dt = (p ■ dp/dt) /. 

Despite numerous independent derivations of equation (1), either on the basis 
of energy-momentum conservation [4, 5] or as the Lorentz force due to the particle’s 
(regularized) self-held [6], it is subject to numerous difficulties; see the recent review 
[7] for an account of these problems and proposed solutions. The most widely used 
alternative to LAD is that introduced by Landau and Lifshitz [8], by treating the self¬ 
force as a small perturbation about the applied force and retaining terms to leading 
order in the small parameter r: 

-a ^ (d,F<^^XbX^ - ^A\F^^F^A . (3) 

m m \ m / 

f The characteristic time r = 2t jZc can be interpreted as the time taken for light to travel across the 
classical radius of the particle, r = q^/dTreowc^. For an electron, Te = 6.3 x 10“^^ s, corresponding 
to Te = 2.8 X 10“^® m. Since radiation damping is proportional to r, radiation reaction effects will 
typically be more prominent for electrons, for which q = —e and m = mg, than for particles with larger 


mass. 
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It is often claimed that (3) is valid provided only that qnantum effects can be ignored, 
and thongh a rigorons demonstration remains elusive there is mounting evidence that 
this is indeed the case [9, 10]. Note that equation (3) can also be presented in terms of 
the electric and magnetic helds, E and B, as [8, 11] 


t = hE+-!^xB 


dt ^ \ 'jm 

dE / p dB \ q 

^7? ^ + — 

dt \7m dt J ym 


E X B + ( — X B 1 X B 

ym 


q. 


y2^2 


(p-E)E+^(p.E)y-^fE+-e-xB') pi, (4) 

\ / y^m^ V ym / | 


ym 


where the total time derivative acting on the E, B helds is d/dt = d/dt + (l/ym)p ■ V. 

Under the conditions expected at ELI, the caveat ‘provided that quantum effects 
can be ignored’ is pertinent. Quantum effects are typically negligible if the electric held 
observed by the particle, E, is much less than the Sauter-Schwinger held [12, 13] typical 
of QED processes, that is provided 


eh !—r-^- E , , 

X := -^V F FacXbX'^ = — -C 1, (5) 

m2 Es 

where Es = mlc^/eh = 1.32 x 10^® V/m is the Sauter-Schwinger critical held. For 
1 GeV electrons in a laser pulse of intensity 10^^ W/cm^ (parameters obtainable at 
ELI), X ~ 0.8 and quantum ehects cannot be ignored. A complete QED treatment of 
radiation reaction is difficult to implement and problematic even to dehne but, provided 
X remains small, a semi-classical modihcation to (3) should be valid [14]. 

An important diherence between the classical and quantum pictures of radiation 
emission can be seen in the radiation spectrum. Classically, a charged particle can 
radiate arbitrarily small amounts of energy at all frequencies. However, in the quantum 
picture, the particle must radiate entire quanta of energy in the form of photons. Thus, 
the energy (frequency) of the emitted photons is limited by the energy of the particle. 
This suppresses emission at high frequencies, and introduces a cut-oh in the spectral 
range of the emitted radiation [14]. As such, it is expected that the ehects of radiation 
reaction are overestimated by classical theories in regimes where quantum ehects become 
important [15], since they consider the particle to be radiating at all frequencies. 

In order to account for this reduction in the ehects of radiation reaction relative to 
the Landau-Lifshitz equation of motion, we follow Kirk, Bell and Arka [16] and scale 
the radiation reaction force by the function g{x)'- 

-a ^ _±pab^^ _ (d,E'^’’XbX'^ - ^A\E^‘^E,dX‘^) . (6) 

m m \ m J 

The full expression for g{x) involves a non-trivial integral over Bessel functions. To 

make this tractable, we use an approximation introduced by Thomas et al. [17], 

^(x) = (l + 12x + 31x'-l-3.7x')""^®. 


( 7 ) 
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It can be clearly seen that, in the classical limit y —)■ 0, we have g (y) —)■ 1, recovering 
the classical eqnation of motion (3). As we move into a more strongly qnantnm regime, 
the qnantnm nonlinearity parameter y increases and the scaling fnnction g{x) decreases, 
in tnrn redncing the effects of radiation reaction. The model essentially rednces to a 
rescaling of the characteristic time of the particle, r —)■ g{x)'r, which can also be applied 
to eqnation (4). For x ~ 1, the stochasticity of qnantnm emission becomes important, 
and the semi-classical model is no longer applicable [18]. At this point, g{x) — 0.18, 
which corresponds to a signihcant redaction in the effects of radiation reaction. 

It is generally accepted that radiation reaction effects will be more readily observed 
in the behavionr of particles than in the radiation they emit [17, 19]. As snch, it is 
important to be able to accnrately determine the distribntion of a bnnch of particles 
evolving according to (3) or its semi-classical extension ( 6 ). Usnally this wonld involve 
solving a Vlasov eqnation [20] or following the evolntion of very large nnmbers of 
particles [ 21 ], either of which is compntationally very intensive. 

In this paper we investigate beam cooling of a particle bnnch dne to classical and 
semi-classical models of radiation reaction. In Section 2 we present a detailed discnssion 
of longitndinal and transverse phase space contraction of the particle distribntion, along 
with an analytical solntion of the classical Vlasov eqnation. The longitndinal particle 
distribntion is introdnced. Since the semi-classical Vlasov eqnation has no analytical 
solntion, in Section 3 we introdnce a new method of accnrately reconstrncting the 
particle distribntion from the trajectories of a relatively small nnmber of particles. 
Classical predictions nsing this method are compared to the analytical solntion with 
excellent agreement. The method is then applied in Section 4 in order to compare 
classical and semi-classical predictions for an electron beam colliding with an intense 
laser pnlse. Finally, we conclnde by snmmarising onr hndings in Section 5. 


2. Particle distribution and phase space contraction 


The evolntion of a particle beam can be described by the Vlasov eqnation for the particle 
distribntion ^{x,u), where = ( 7 , u) is the 4-velocity. Position and velocity are 
considered as independent phase space variables. The Vlasov eqnation for ^ can be 
expressed as 



d^ 

ds 



V 


0 , 


( 8 ) 


where V is the phase-space volnme element and (ds{.x, u) describes the rate of change (he. 
expansion or contraction) of V with proper time s. (Technically, fdg is the phase-space 
divergence of the vector held X = u°‘d/dx°‘ + Afdjdu^ (where A is the acceleration) 
associated with the how d/ds, given by the Lie derivative LxL = /d^V, see [20].) Capital 
Latin indices take three valnes. Unlike the Lionville eqnation (or the case with no 
radiation reaction) the phase-space volnme element is not preserved by the how, fig 7 ^ 0 . 

To facilitate investigation of the interaction of a particle bnnch with a laser pnlse. 
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we introduce the (null) wavevector k such that the phase of the pulse is 


(f = —k ■ a; = cut — k • X. (9) 

The orthogonal (transverse) vectors e, A satisfying 

= 1 and /c-e = /c-A = e- A = 0, (10) 

together with k and the null vector £ (defined to satisfy £ • e = £ • A = 0 and k ■ I = — 1) 
form a basis. In addition, the coordinates 

f = e ■ X, a = X ■ X and = —i ■ x (11) 


are also defined, along with the corresponding velocities u^j), and m,/,. However, 

is not independent and may be found from the normalisation condition u°‘Ua = 
— 2u^u^ = — 1. We note that Greek subscripts are used only as labels and are 
not free indices. 

For a plane wave with arbitrary polarisation, the electromagnetic field tensor F 
depends on spacetime only through the phase 0 , and takes the form 

= a,(0) {e^h - k^e,) + ax{(t>) {X^h - k^Xf ), (12) 

Tflj 

where the functions are dimensionless measures of the electric field strength in the 

e,A direction. The corresponding electric and magnetic fields are E = {mu/q)[ae{(j))e + 
aA(0)A] and B = kx E/cn, where the orthogonal unit 3-vectors e, A satisfy k-e = k-A = 

0 . 

In a similar manner, we assume that the particle distribution also depends on 
spacetime only through the phase 0, such that ^{x,u) = ^(0, Mg-). The Vlasov 

equation is then written 


Vj(j) 




A 




jds 


(13) 


where G {m^, Mg, Wo-}, and the accelerations A^ G {Al^, Alg, A^j} follow from the single¬ 
particle equations of motion. Dividing through by Mg,, we have 

(9 fA^\ 

whe.e fS-^.{-). (14) 

The quantity fd is responsible for any phase space contraction {fd < 0) or expansion 
{(d > 0 ) of the particle distribution, and the associated change in electron entropy [ 22 ]. 

For a highly relativistic particle beam colliding with a laser pulse (the scenario in 
which radiation reaction effects are most prominent), we are mainly interested in the 
dependence of ^ on 0 and Mg,. An advantage of the coordinate system (9)-(ll) is 
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that it decouples the longitudinal from the transverse velocity in the Lorentz invariant 
measure, d^x/'j = du^du^du^/u^. Hence we can dehne the longitudinal distribution 


= \ ^ du^dua, 


(15) 


which satishes the reduced Vlasov eguation 
^+Al/ = 0, where 


ou^ \u^ 


(16) 


Here, /3\\ describes the longitudinal phase space contraction. The transverse contribution 
is then 


(3^ = (3- /?|| 


du^ \u^) dua \u^ J ' 


(17) 


Note that this reduction to the longitudinal distribution is purely a consequence of the 
coordinate system, and does not rely on the plane wave assumption. 

It is at this point that a decision must be made as to the appropriate single-particle 
equations of motion. While there are many classical models for radiation reaction [7], we 
start by considering the Landau-Lifshitz equation given by (3), before moving on to the 
semi-classical extension (6). This is in part motivated by the existence of an analytical 
solution to the single-particle Landau-Lifshitz equation [23]. In our coordinates, the 
Landau-Lifshitz equations in the plane wave (12) are 


= -rgful 

= —u^{a^ + Tu^a'^ — (18) 

A„ = -u^{ax + Tu^a'x) - rafulu^, 

where gf{(f)) = and prime denotes differentiation with respect to (j). 

Inserting these equations into (16) and (17), we hnd for the classical case 

/3|| = = -2ra\^ < 0. (19) 


It is immediately apparent that half the contraction of the distribution occurs in the 
longitudinal and half in the transverse directions. 

The semi-classical equations of motion are just (18) with the replacement r —)■ 
g{x)'T- However, since = 3Ta{(f))u^/2a (where a is the hne structure 

constant) depends on (but not on the transverse velocities) we pick up an additional 
contribution to the longitudinal phase space contraction: 


a 


9^ + 


dg A^ 

dU(j) U(i) 


and 


Ai = g^w + 


dg A<i>_ 

dUfj) Ufj) 


( 20 ) 


whereas, the transverse contraction is simply scaled by g{x)'- 

P±= !}- P\\= 30 - Al) = Sffx = sX 


( 21 ) 
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X X 


Figure 1. Reduction of longitudinal beam cooling. Part (a): the ratio (23) of the 
longitudinal to the transverse phase space contraction in the semi-classical model. The 
dashed line shows the classical ratio /3||//3 _l = 1- Part (b): ratio of the semi-classical 
longitudinal contraction to the classical Landau-Lifshitz result. 


Thus, as quantum effects become more important and g{x) decreases, the semi- 
classical model predicts a reduction in both the longitudinal and transverse phase space 
contraction (reduced beam cooling). As well as this scaling of the classical contraction 
by g{x)y there is an additional longitudinal heating given by 

dg _ dg dx 
dU(p U(p dx ldU(p U(i) 


gX//"(x)(l2 + 62x+ 11.1x2) 


> 0 , 


( 22 ) 


where (3± = gj3± = —2ra2(0)5((x)u</, = —4Q;a(0)x5'(x)/3- The ratio 

/3± 2d\ogx 

= 1 - ‘^xg^^^x) (12 + 62x + ll.lx') < 1 (23) 

measures the strength of the longitudinal compared to the transverse phase space 
contraction. This is shown in Fig. 1(a) for the interval x ^ [0, !]• Even for the weakly 
quantum regime in which the semi-classical model remains valid, we observe a signihcant 
reduction in longitudinal beam cooling. This is especially clear when comparison is made 
with the classical result /3\\ as shown in Fig. 1(b). We see that where x = 0.2 there 
is nearly a 60% reduction in the longitudinal contraction experienced compared to the 
Landau-Lifshitz model. 

For the case of the classical Landau-Lifshitz theory in a plane wave, the Vlasov 
equation (14) may be solved analytically for the particle distribution: 


^(0, u^, u^, u^) = ^ {ctP, mJ, ul) (24) 

where {0°, } are the initial phase and velocities of a particle with 

at phase 0. In a similar manner, the longitudinal distribution is found to be 




(25) 
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The contraction/expansion of phase space is contained in the function 


= T dd 


(26) 


Solutions to equation (18) [23] can then be used to rewrite {0°, } in terms of 

the independent variables {0 , Mo-}: 


,,o _ 


1 - tu^G{4>) 

- TU^(yVe{(p) + CteiG) - ae(0°) 






1 - tu^G{4>) 

Ua - TuJWx{(t)) + aA(0) - ax{(tP) 


where the functions 


1 - TU^GiG) 


G{G)= / dd gf{d) 


Vi(0) = dd ai{d) 


= / dd a,{d)G{d) 


H( 0 ) 

-Va(0), 


(27) 


(28) 


with i G {e, A} depend only on the properties of the laser pulse. 

Using equations (27) we can express u^{d) in equation (26) in terms of the 
independent variable 


A(()),m^) = TU^ / dd 


gf{d) 


1^0 I-Tu^{G{(p) - G{d)) 

= - In - TU^GiG)^, 
such that the distribution becomes 


^ (0, W(t) 


(1 - Tu^GiG)) 


4 ’ 


(29) 


(30) 


or for the longitudinal distribution 


/ ( 0 ) 


/ l-r^Gw) 


[1-tu^G{G)) {^-ru^G{G)y 


(31) 


This latter result is in agreement with observations made by Neitz and Di Piazza [24], 
and we see that the longitudinal distribution is only sensitive to the properties of the 
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laser pulse through the function After the pulse has passed, Q becomes constant 

and is proportional to the fluence of the pulse. Final-state properties of the longitudinal 
distribution therefore depend only on the total energy contained in the laser pulse, and 
are insensitive to how that energy is distributed within the pulse. The full distribution 
on the other hand, depends additionally on the integrals Wj given in equation (28). 

Although the reduced Vlasov solution (31) is somewhat simpler than (30), and 
captures the key features of the electron beam itself, the solution is not sufficient to 
calculate the transverse current density, and hence cannot be coupled to Maxwell’s 
equations to determine the radiation produced by the electron beam. However, if 
the transverse momentum spread is sufficiently small, we can approximate the full 
distribution by 


^(0, U^, U^, Ua) = /(0, U^)6 - A;(0, u^)) 6 - A’a(0, u^)) , (32) 


where the ^-functions restrict the transverse velocities to the submanifold A). Then, in 
addition to (16), equation (14) yields 


d(j) 


- rgful 


du^ 


{di -\- TU(j,al + TU^^Xi) , for i G {e. A}. 


(33) 


Note that (33) indicate that the distribution is concentrated on a surface in phase space 
that itself satishes the Landau-Lifshitz equation. 

Given solutions to the reduced Vlasov equation (31) and the transverse Landau- 
Lifshitz equation (33), the current can be written 


r = q [ du^du^^ = Jl + Qt -f Jl (34) 

J ^0 

with jf_ and g evaluated as 

Jl = f f + qX'^ [ f Xx— and g = q [ fdu^ . (35) 

J J J 

We could also calculate directly, but it follows more straightforwardly from charge 
conservation, daj^ = 0 . 

In the following, we restrict our attention to the longitudinal distribution 
and longitudinal beam cooling, as this is more readily measurable in experiments than 
the transverse cooling. However, the transverse cooling, which can be considerably 
greater, can be determined from equation (23). 


3. Numerical (re)construction of the particle distribution 

The motion of a single charged particle colliding with a laser pulse, including radiation 
reaction, has been extensively studied [9, 23, 25, 26]. As shown in Section 2, the Vlasov 
equation with classical readiation reaction can be solved analytically. However, this 
is not the case for the semi-classical model ( 6 ) or for stochastic models of radiation 
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reaction in the quantum regime. Instead of attempting to solve a Vlasov-type equation 
on the phase space numerically, which would require significant computing resources, we 
propose an innovative method which allows for the dynamics of a particle distribution 
to be explored using single-particle equations of motion such that the distribution can 
be efficiently reconstructed. While this approach is quite general and could be used for 
a variety of systems, here we consider a distribution of particles subject to equation (3) 
and its semi-classical extension (6), without particle-particle interactions§. 

Assuming that the laser pulse can be approximated by a plane wave with compact 
longitudinal support ||, any spatial spread in the initial particle distribution would only 
determine the moment when each particular particle enters the pulse. For simplicity, 
we therefore take all particles to originate from the same point. This is reasonable as 
we are primarily interested in the longitudinal momentum distribution. 

Since our pulse is modelled by a plane wave and we focus on the longitudinal 
properties of the distribution, we consider the initial momenta to be strongly peaked 
about zero in the transverse directions. As such, the initial distribution can be taken to 
be a Maxwellian distribution for the (longitudinal) momentum p (in units of me) 


/(0 = 0 , p ) 


Np 

, - exp 


{p - p)^ 

29 


(36) 


with 0 = cut — k ■ X the phase, 6 the variance of the distribution, and Np the number of 
particles. The momentum p is related to our velocities of Section 2 by 


p = 


where 


l + u^^ + ul + (u^/uj)^ 
2{u^/u) 


(37) 


We stress that this initial distribution is chosen for its simplicity; alternative 
distributions could be used where appropriate (such as Maxwell-Jiittner). 

Typically, one would sample the distribution at random, which would require a 
large number of particles to accurately represent the distribution. Instead, since the 
particle number is simply 

/ OO 

dpfi<P,p), (38) 

■OO 

we determine the momentum spacing dp between the particles from the initial 
distribution by truncating the integral in (38) so that the particle number increases 
by unity in the given momentum interval: 


1 


^ j dp f(0,p) ~ f(0,p)Sp. 


(39) 


§ For a highly relativistic particle bunch, these interactions can be neglected on the time scale of the 
laser interaction. 

II A function has compact support if it is zero outside a finite interval. 
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This leads to a set of Np = 2Nc + 1 initial momenta 1^(0) = {pi(0)} for i e [—Nc, A^c], 
with the Pi generated iteratively from po = P and p±i = p ± l//(0,p) using 


Pi = Pi- 2 i + 777^ - r with ^ = sgn(i). (40) 

/ lO, Pi— 

The momentum space is not sampled uniformly, instead more particles are located in 
regions where the distribution function is large. 

As the evolution proceeds, this procedure is applied in reverse to reconstruct the 
distribution. The set of momenta V (0) is ordered such that pi+i > pi and used to find 
Spi{(j)) = (pj+i(0) — Pi_i(0))/2. The velocity distribution is then defined to be 

' Spi{(j)) 

Reconstruction of a distribution from a particle sample can be problematic, but in our 
formalism it becomes quite natural. This is achieved by using the momentum spacing 
between particles to determine the value of the distribution such that equation (39) is 
satisfied for all 0 (he. integration over each of the measured momentum spacings always 
contributes a single particle to the total particle number). The closer the measured 
momenta are together, the ‘more likely’ one is to have a particle in that momentum 
range, resulting in a larger value for the distribution. 

The dehnition (41) allows properties of the distribution to be calculated directly 
from the momenta of the individual particles. The mean is simply evaluated as 



(41) 


p(0) = (p(0)) = ^ 


dpp/(0,p) 


iVp ^^ 

* =1 


(42) 


In a similar manner, higher-order moments of the distribution may be calculated 
straightforwardly as: 

^n(0) = ([p-M0)]”> - [p*(0) (43) 


For example, the variance 6*(0) = ^ 2 ( 0 ) and skewness 5'(0) = A'3(0)/X2'^^(0). We note 
that angular brackets (■ ■ ■) will be used to denote an average over the bunch, not time. 

A major advantage of this new methodology is that the distribution can be 
accurately represented and reconstructed using far fewer particles than would be 
required using random sampling. Figure 2(a) shows a simple Gaussian distribution 
with zero mean and unit variance, A{z) = (1 /a/^) exp(— 2 :^/ 2 ), reconstructed using two 
methods. First, the distribution A{z) was sampled according to the iterative relation 
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Figure 2. (Colour online) Part (a): A Gaussian distribution A{z) is reconstructed 
using the new method described by equations (40) and (41) with Np = 401 particles 
(red, solid) and compared to the standard approach of random sampling, using both 
iVj = 401 (blue, dot-dash) and = 100,000 (black, dashed) particles. Part (b): 
Variation of the measurement of the initial relative momentum spread di as the number 
of particles Np is increased, compared to the desired value, dd- 


(40) with Np = 401 and then reconstructed using (41) (—). Next, = 100,000 
random numbers were generated from A{z) using a Mersenne-Twister pseudo-random 
number generator, from which the distribution was found using 100 fixed-width bins 

(-). Even with 100,000 random samples, the standard approach does not describe 

the Gaussian perfectly, while the new method presented in this paper does an excellent 
job using only 401 points. For comparison, a sample of = 401 was also reconstructed 
using 50 fixed-width bins (-■-). The advantage of the new method with such a small 
sample size is clear. 

Figure 2(a) also shows how sampling the distribution according to equation (40) 
cuts off the tails of the distribution, which will have an effect on the measured properties 
of the reconstructed distribution. The finite number of particles used to represent the 
distribution causes the measured relative momentum spread (defined by equation (44) 
below) to be less than that specified when defining the initial distribution. Essentially, 
it comes down to the in equation (39), compared to the definition given by equation 

(41) . As Np is increased, the approximation in equation (39) improves, and the 
measured value approaches the desired value. In addition, with more particles the 
distribution is sampled further into the tails. Figure 2(b) shows the measured initial 
spread, d,, as a fraction of the desired spread, dd, when the particle number is varied 
from as few as 11 up to Np = 2001. We see that the approximation quickly improves 
as Np is increased up to about 500. The value Np = 401 is chosen to give less than 
0.5% error in the initial measured momentum spread. In practice, good agreement can 
be found with lower Np, with the caveat that properties sensitive to the tails of the 
distribution (such as the skewness) may be strongly affected. (This is conhrmed in Fig. 
3, where the distribution and its statistics calculated from the analytical solution to the 
classical Vlasov equation are compared to those obtained with this new method.) 
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4. Interaction of a particle bunch with high-fluence laser pulses 


The analytical solution to the Vlasov equation including radiation reaction according 
to the Landau-Lifshitz theory derived in Section 2 predicts that the collision of an 
energetic electron beam with a high-intensity laser pulse leads to a signihcant contraction 
of the particle phase space, resulting in a reduction in the relative momentum spread. 
Agreement between this analytical solution and numerical results obtained using the 
approach discussed above is shown below to be excellent. As previously observed, 
classical beam cooling depends only on the total fluence of the pulse, rather than 
its duration or peak intensity independently [24, 27]. However, for the semi-classical 
extension, the Vlasov equation is no longer tractable. This highlights the value of 
our approach and, to demonstrate the use of our proposed method in such a case, we 
consider the importance of quantum effects in the interaction of an electron bunch with 
a high-intensity laser pulse. 

To establish the impact of quantum effects on the evolution of the particle 
distribution subject to radiation reaction, we introduce the relative momentum spread 
and the momentum skewness (calculated from the mean p and variance 6): 


(T(0) 


VW) 

p{(j)) 


and S'(0) 


03 / 2 ( 0 ) 


(44) 


The former gives a measure of the beam quality, while the latter indicates how symmetric 
the distribution is about its mean. 

We restrict our attention to a linearly polarised iV-cycle plane wave pulse (12), 
modulated by a sin^-envelope [9], with = 0 and = a(0), where 


a(0) 


Qq sin((;/)) sin^ {ircj)/L) for 0 < 0 < L 
0 otherwise 


(45) 


where Oq is the dimensionless (peak) intensity parameter (the so-called “normalised 
vector potential”) and L = 2ttN is the pulse length^. This pulse shape offers compact 
support, allowing the particles to begin and end in vacuum. The total ffuence (energy 
per unit area) of the pulse is proportional to 

L= df) a‘^{(j)) = — Nal- (46) 

Jo 8 

In this work, S is kept constant, which fixes oq for each N. It has been shown [24, 27] 
(see also Section 2) that the classical Landau-Lifshitz prediction for the final state of a 
particle distribution emerging from the pulse is completely determined by the ffuence, 
whereas quantum effects are expected to depend on the value of oq itself. We are then 
able to explore the impact of the reduced emission in the quantum model with varying 

^ For the sin^-envelope used in this work, the full-width half-maximum (FWHM) duration is half of 
this value. 
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Figure 3. (Colour online) The phase space evolution of the distribution function 
/(</), p) predicted by the analytical solution (31) of the reduced Vlasov equation (16), 
compared to numerical results obtained using the new method presented in Section 3. 
Results are presented for (a) N = 20 and (b) iV = 5 cycles, where L — 2'kN is the 
pulse length. The fluence has been kept constant, with iVog = 9.248 x 10^. 

Oq while maintaining the same classical prediction. This allows us to explore the relative 
importance of quantum effects. 

To motivate this study, parameters have been chosen to be relevant at the 
forthcoming ELI facility. We have chosen to consider Noq = 9.248 x 10^ which, for 
= 20 with a wavelength of A = 800 nm, represents a full-width half-maximum pulse 
duration of 27 fs with peak^ intensity 2 x 10^^ W/cm^. We have investigated pulses of 
length N G [5, 200] cycles (together with their corresponding uq) counter-propagating 
relative to a bunch of Np = 401 particles, with an initial momentum spread of 20% 
around a/ 1 -b = 2 x 10^. This corresponds to an average particle energy of just 
over 1 GeV, which should be well within the capabilities of the laser-plasma wakeheld 
accelerator at ELI. 

Before comparing predictions of the classical and semi-classical models, we briefly 
conhrm the validity of our method by comparing numerical results with the analytical 
solution (31) obtained in Section 2. Figure 3 shows the interaction of a 1 GeV electron 
beam with a plane-wave laser for two pulse lengths, iV = 20 (uq — 22) in part (a) and 

+ Peak intensity is obtained from /peak = (47r^m.gC®eo/e^A^)aQ ~ 2.74 x 10^°(ao/A)^ W/m^. 
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N = 5 (ao ~ 43) in part (b). The left-hand panels show the nnmerical resnlts obtained 
using the new method described in Section 3, while the right-hand panels show the 
solution (31) using the initial distribution f(0,u^) corresponding to f(0,p) given by 
(36) with p = \{Ufj,/uj — u/u^). The momentum p is evaluated during the evolution 
using (37), along with the solutions (27) for when = 0. The agreement 

is excellent. Note that the measured values for the initial and final momentum spread 
also agree, while the skewness is underestimated by the numerical method (as discussed 
in Section 3). 

Figure 4 shows the variation of the particle distribution on the (0,p) phase space. 
As can clearly be seen in moving from the classical Landau-Lifshitz theory (left) to 
the quantum model (right), there are noticeable differences in the mean p, spread a, 
and skewness S of the distribution. We first note that the deficit in measuring the 
initial dj = 19.9% < 20% is due to the finite number of particles used to represent the 
distribution (as discussed at the end of Section 3 and illustrated in Fig. 2(b)). 

For the classical theory, the hnal distribution only depends on the fluence of the 
pulse, though this does not prevent the system from taking different routes along the 
way. As the number of cycles is decreased, very different intermediate behaviour is 
observed in Fig. 4, yet the measured properties of the final distribution support this 
prediction: in each case, we measure the mean momentum pf = 1197.7 with a relative 
spread d/ = 12.5%. This represents a significant contraction of the phase space, where 
the average energy of the particle bunch decreases significantly, as does its thermal 
spread (beam cooling), and the distribution becomes more sharply peaked. In addition, 
we find the development of a negatively-skewed distribution with Sf = —0.46. In the 
classical model this is readily understood, since the higher a particle’s momentum the 
more it radiates. This causes particles in the positive tail of the distribution to be slowed 
down more than those in the negative tail. 

The introduction of a semi-classical model in which the effect of radiation reaction is 
reduced by the function g{x) given by equation (7) results in a reduction in the amount 
of phase space contraction. Figure 4(a) for N = 200 clearly demonstrates this, with 
the hnal average momentum pf = 1301.1 only slightly higher than the classical case. 
The hnal relative momentum spread is now 14.4%, showing that the hnal distribution 
is less sharply peaked. While remaining negative, the skewness reduces in magnitude to 
—0.280, because it is precisely the higher-energy particles (which were classically most 
ahected by radiation reaction) that now have this damping suppressed due to larger y 
(smaller g{x))- 

These changes become more pronounced as we move to higher intensities (by 
reducing the number of cycles). For N = 20, as shown in Fig. 4(b), we hnd that 
Pf = 1451.8 and dj = 16.6% have both increased, with the skewness also increasing to 
Sf = —0.129. This trend continues to A^ = 5 as shown in Fig. 4(c). In this case, very 
little beam cooling occurs for the quantum model, with the final relative momentum 
spread taking the value dj = 18.0% around pf = 1581.3. The profile also remains much 
more Gaussian, with Sf = —0.0597. 
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Figure 4. (Colour online) The phase space evolution of the distribution function 
/(</), p). Classical predictions are shown in the left-hand panels, while the corresponding 
semi-classical results are presented to the right. Values of the initial and final relative 
momentum spread and momentum skewness are displayed in each figure. The pulse 
length is reduced from N = 200 cycles in part (a) to iV = 20 cycles in part (b), and 
finally to V = 5 cycles in part (c). In each case, we observe an increase in the final 
mean momentum and its spread (reduction in beam cooling) predicted by the quantum 
model. 


The reduction in phase space contraction (beam cooling) observed here is in 
agreement with previous predictions [28], which provides further validation of our 
method for reconstructing the particle distribution. Using this new method, it has 
been possible to investigate the effects of the semi-classical model on a distribution of 
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Figure 5. (Colour online) Part (a): variation of the final relative momentum spread 
difference Saf as a percentage of the total classical change, The evolution of the 

average quantum parameter (x) according to the classical and semi-classical theories 
are plotted in part (b) for N = 20. 


particles. The distributions are nicely reconstructed and do not feature any artefacts, 
in contrast to other approaches [17], which emphasises the power of our method. 

Figure 4 also shows how the difference between the classical and quantum results 
increases as the intensity is increased (or as N is decreased). It is therefore interesting 
to consider the difference 6af = aj^ — af as a fraction of the total (constant) classical 
change in momentum spread, = dj — d^h This can be found in Fig. 5(a), where we 
see that for N = 5 the two predictions differ by about 75%. As N is increased, this ratio 
is reduced because the average quantum parameter (y) becomes smaller and radiation 
reaction is not so heavily suppressed. It would be expected that the two models converge 
as A —)■ oo. 

In cases where there is a large discrepancy between the predictions of the two 
theories, it is especially important to be confident in the validity of the model. As a semi- 
classical model, we expect it to remain valid into the weakly quantum regime, such that 
particles experience instantaneous values -C 1. Figure 5(b) shows the evolution of 
the bunch-average (y) as the bunch moves through a laser pulse with A = 20 according 
to the classical and semi-classical models. Initially, there is good agreement between 
the two models, until the bunch approaches the centre of the pulse, where the intensity 
becomes higher and the models are signihcantly different. For completeness, we note 
that our highest intensity case with A = 5 satished (y)^ < 0.22. 


5. Conclusions 

The next few years will see the emergence of a number of new high-power laser facilities 
operating at unprecedented held strengths, providing access to fundamentally new 
physical regimes. This will allow us to experimentally probe previously untested areas 
of physics, such as the long-standing question of radiation reaction. 

In this paper, we have analysed the transverse and longitudinal cooling of a 
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relativistic electron beam as it interacts with an intense laser pulse, according to classical 
and semi-classical theories of radiation reaction. In the classical theory, we have found 
these two contributions to be equal, but quantum effects break this symmetry, leading 
to significantly less cooling in the longitudinal than the transverse directions. 

To facilitate evaluation of the longitudinal beam cooling effects, we have introduced 
an innovative method to efficiently and accurately calculate the distribution function 
for an electron beam interacting with an intense laser pulse. This has been validated 
by comparison with an analytical solution to the Vlasov equation in the classical case, 
and used to compare classical and quantum predictions of radiative cooling. We have 
found that quantum effects can significantly alter the beam properties and, unlike the 
classical case, can be influenced by the shape of the laser pulse, not just its energy. 

As we move into the quantum regime where final-state electron beam properties 
become sensitive to pulse shape, it is becoming increasingly important to have an 
efficient method in order to investigate the full parameter space. The approach 
developed here to facilitate this study of beam dynamics provides a powerful tool with 
wide-ranging application within the discipline. 

The results presented in this paper are limited to the semi-classical case -C 1. 
However, it should be noted that, for the longitudinal beam cooling, this restriction is 
due to the use of a deterministic equation of motion, and not the method of sampling 
and reconstructing the distribution. There should be no obstruction to exploring more 
strongly quantum regimes (such as higher initial beam energies ~ 5 GeV available at 
ELI) using this approach with a stochastic equation where photon emission probabilities 
are determined by strong field QED, as in [29, 30]. This will be addressed in future work, 
along with an investigation of stochastic transverse beam cooling. 
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